Terahertz tunable three-dimensional photonic jets

Highly localized electromagnetic field distributions near the “shadow-side” surface of certain transparent mesoscale bodies illuminated by light waves are called photonic jets. We demonstrated formation of three-dimensional (3D) tunable photonic jets in terahertz regime (terajets, TJs) by dielectric micro-objects -including spheres, cylinders, and cubes-coated with a bulk Dirac semimetal (BDS) layer, under uniform beam illumination. The optical characteristics of the produced TJs can be modulated dynamically through tuning the BDS layer’s index of refraction via changing its Fermi energy. It is demonstrated that the Fermi energy of BDS layer has a significant impact on tuning the optical characteristics of the produced photonic jets for both TE and TM polarizations. A notable polarization dependency of the characteristics of the TJs was also observed. The impact of obliquity of the incident beam was studied as well and it was demonstrated that electromagnetic field distributions corresponding to asymmetric photonic jets can be formed in which the intensity at the focal region is preserved in a wide angular range which could find potential application in scanning devices. It was found that the maximum intensity of the TJ occurs at a non-trivial morphology-dependent source-angle.


Geometry of the simulation
We performed finite-difference time-domain (FDTD) simulations in 3D using Lumerical software package (Lumerical Inc., Canada). Figure 1 shows the schematic of the proposed structures, a sphere, cylinder, and cube made of high-density polyethylene (HDPE) with refractive index n 1 = 1.51 at the THz frequency range 80 , surrounded by air (n bg = 1) and covered with a uniform BDS layer of thickness d DS and index n DS .
The quantitative key parameters of a typical PNJ or TJ are schematically presented in Fig. 1.Focal distance (FD) was defined as the distance from the proximal surface of the generating object to the focal point (locus of the maximum intensity) of the TJ.The TJ's transverse size, quantified as the full-width at half-maximum (FWHM) in the xy-plane along y-axis, was calculated at the focal point.The TJ's length was defined as the distance between proximal and distal 1/e contour of the intensity along the beam direction (x-axis).
The dielectric function, ε DS , of the BDS can be written as 67 : in which f is the frequency of the incident wave, ε b is the background's effective dielectric constant, ε 0 is the permittivity of vacuum, i = √ −1 , and σ DS f , E f represents the BDS's conductivity at a given frequency and Fermi Energy ( E f ).The following equations give the complex optical conductivity of the BDS 67 : in which k F = E F / v F is the Fermi momentum, v F = 10 6 m/s is the Fermi velocity, ε = E/E F , ε c = E c /E F , E c is the cut off energy above which the dispersion relation of the BDS becomes non-linear ( ε c = 3 is our case), and T = 300 K is the temperature.Symbols k B , e, , and g represent the Boltzmann constant, an electron's charge, the reduced Planck constant, and the degeneracy factor (g = 4), respectively.
A uniform beam propagating along the x-axis, with unit amplitude, whose electric field oscillates along the i) z-axis (TE) and ii) y-axis (TM) was incident upon the proposed structures to investigate the dependency of the produced TJs on the polarization of the light source.The computational mesh size was Δx = Δy = Δz = 0.3 μm (~ λ/140-λ/100) to secure converging results.The perfectly matched layers (PMLs) boundary condition (34-μm-thick) were employed around the simulation region to prevent spurious reflection 81 .The incident uniform beam had a finite size encompassing the whole cross-sectional area of the object; its size along y-axis (z-axis) is 1200.3μm (2201.4μm) for sphere, 694.4 μm (680.9 μm) for cylinder, and 481.3 μm (691.8 μm) for the cube.The distance between the center of the source and center of the all dielectric sphere, cylinder, and the cube objects is 30 μm.Each run took about 8-9 h, 4-5 h, and 4-5 h for a sphere, cylinder, and cube, respectively, in our system with 64 GB RAM.
We numerically investigate the ability of the proposed BDS-coated dielectric objects (sphere, cubic and cylinder) to produce TJs, and control their optical properties by altering the BDS layer's Fermi energy at room temperature.The BDS layer is considered to have similar properties to Cd 3 As 2 (cadmium arsenide) and Na 3 Bi (trisodium bismuthide) with ε b = 12 and g = 4 67 .In contrast to 2D Dirac fermions in graphene or on the surface of 3D topological insulators, topological DSs possess 3D Dirac fermions in the bulk.Cd 3 As 2 is a 3D DS material with separated Dirac points and linear energy dispersion in momentum space.It is a tetragonal material with space group I41/acd and lattice constants of a = b = 1.26 nm and c ≈ 2.54 nm, which remains nearly cubic ( 2a = 2b ≈ c ) 82 .This cell contains 96 cadmium atoms and 64 arsenic atoms with two cadmium vacancies.The Dirac equation describes the massless electrons in Cd 3 As 2 , suggesting a conical band that has twice degenerated due to spin.The crystal structure of Na 3 Bi, as a topological DS, is a body-centered tetragonal with space group of P4/mmm.The electronic band structure measurements and calculations of Na 3 Bi DS show a cone-shape with linear dispersion 68 .The lattice constants of Na 3 Bi at the ground state are a = 3.4116A • and c = 4.9530A • 83 .
Cd 3 As 2 , as a 3D BDS, has attracted research interests due to its chemical stability and extraordinary optical response 84 .It has been fabricated in thin films and at nanoscale through different techniques, including physical vapor deposition (PVD) 85 , pulse laser deposition (PLD) 86 , molecular beam epitaxy (MBE) 87 , chemical vapor deposition (CVD) 88 , and self-selecting vapor growth (SSVG) 82 .
The complex index of refraction of a non-magnetic medium is obtained through ñ(ω) = √ ε DS (ω) .The behavior of the real and imaginary index of refraction of the considered BDS versus the Fermi energy ( E F = 5-60 meV) are plotted in Fig. 2 at different frequencies (6-10 THz).It can be seen that as the frequency decreases, the amount of change in the refractive index (dynamic range) over the range of Fermi energies increases, indicating higher degree of tunability of the produced TJs.However, at 6 THz, the real part of the index drops quickly to low values close to and smaller than unity.Hence, f = 7 THz (λ = 42.83μm) which maintains a reasonably high index while providing the highest dynamic range was selected for all simulations in this work.The real part (n r ) of the index decreases gradually with Fermi energy up to a certain energy beyond which it drops suddenly.The imaginary part (n i ) of the refractive index decrease as the Fermi energy increases up to a certain energy beyond which it starts to increase.We chose E F = 10 meV as the lower end at which n r = 2.669 and n i = 0.071.At 7 THz and for E F > 55 meV, the real part of the index reaches unity, so E F = 55 meV was selected as the higher end at which n r = 1.187 and n i = 0.066.Due to a relatively low n i , the absorption associated with the BDS layer will be low in the selected Fermi energy range ( E F ∈ 10-55 ).
It should be mentioned that in practice the E F and consequently the BDS's index can be modulated dynamically through applying a voltage [89][90][91][92] , or by an alkaline surface doping 68,93 .In this way, the propagation of light can be modulated by the BDS-coated objects, affecting the generated TJs' characteristics.
In order to demonstrate the formation of tunable 3D TJs, we considered three cases: (i) a BDS-coated dielectric sphere with radius R = 500 μm (~ 11.7λ) and thickness of the BDS layer d DS = 1.5 μm (~ 0.035λ); (ii) a cylinder with radius R = 300 μm (~ 7λ), height h = 180 μm (~ 4.2λ), and thickness of the BDS layer d DS = 1.5 μm; (iii) a cube with side length R = 300 μm and thickness of the BDS layer d DS = 1.5 μm.The thickness of the BDS layer was kept small in all cases to minimize the absorption.The object's size for each case was obtained through simulating objects with sizes in the range 2λ up to 25λ, where λ is the wavelength, with 0.1λ increments and selecting the size that produced the highest intensity TJ.This is the size range over which production of photonic jets is more efficient 53 .Next, the BDS layer with different thicknesses between 0.5-2 μm, with 0.1 μm increments, was added and the optimum thickness was selected.Here, there is a trade-off between refractive power of the layer and its absorption imposed by the imaginary component of the refractive index.Next, we studied the impact of the obliquity of the incident uniform beam on the TJ formation in the abovementioned cases.Figure 1b,c show the geometry of the oblique incident case on a sphere/cylinder and a cube, respectively.The input source angle was defined as the rotation angle of the source along x-axis as represented by θ in in Fig. 1b,c.The emergence angle of the TJ was defined as the angle between the x-axis and the line passing through the location of the maximum intensity outside of the object and the center of the object as shown by θ out in Fig. 1b,c.Other parameters of the simulation remained unchanged.

Results and discussions
The distribution of the electric field of the TJs produced by the BDS-coated sphere at Fermi energies of E F = 10 meV and E F = 55 meV, for TE-polarized incident wave, are presented in Fig. 3a,b, respectively.The corresponding results for the TM-polarized case are presented in Fig. 3c,d, respectively.The modulation of the characteristics of the produced TJs, realized through changing the BDS layer's Fermi energy, are presented in Fig. 3e,f.Since the intensity is proportional to the second power of the electric field distribution, in order to calculate the TJ's FWHM, the location of the maximum intensity is identified first.Next the intensity profile along a line parallel to the y-axis passing through the location of the maximum intensity is obtained.We observed a decreasing (increasing) trend in FD and length (intensity and FWHM) of the TJs as the Fermi energy increases.For example, for the TM (TE) incident wave, the FD decreased from 68.4 (74.7) μm to 62.8 (63.2) μm when the E F of the BDS layer was increased from 10 to 55 meV, providing 5.6 (11.5) μm dynamic range.Also, within the same Fermi energy range, the enhancement in intensity of the TJ was increased from 313.5 (332.5) to 572.9 (624.4) for the TM (TE) case.The FWHM decreased from 35.4 (32.6) μm to 31.23 (26.4) μm while the TJ length increased from 78.4 (84.6) μm to 88.9 (97.2) μm in which the corresponding amounts of tuning were 4.17 (6.2) μm and 10.5 (12.6) μm, respectively.It can be seen that the TE case provides more dynamic range for the tuning compared to the TM case.The focal distance and TJs are longer, more intense, with smaller FWHM in TE cases compared to the TM cases.It has been suggested that TM polarization generates wider PNJs compared to the TE polarization due to longitudinal component of the TM-polarized incident light 94 .
The distribution of the electric field of the TJs generated by the cylinder at E F = 10 meV and E F = 55 meV, for TE polarized incident wave, are shown in Fig. 4a,b, respectively.The corresponding cases for the TM-polarized case are presented in Fig. 4c,d, respectively.The plots of the TJ's FD and intensity, as a function of E F , are presented in Fig. 4).The plots of the TJ's FWHM and length, as a function of E F , are presented in Fig. 4f.It can be seen that the optical characteristics of the produced TJs follow a similar trend as in the corresponding sphere cases in terms of the Fermi energy.However, the sphere provided TJs with an order of magnitude more intensity.The intensity enhancement of the TJ varied between 28.9 (27.Using geometrical optics paraxial approximation 18 , we can calculate the focal distance of a core-shell system with circular cross-section with the core and shell radii of R c and R s and refractive indices of n c and n s , respectively , as: www.nature.com/scientificreports/For our core-shell spheres with R c = 500 μm and R s = 501.5 μm, we have FD = 242.8μm and 239.7 μm at E F = 10 meV and E F = 55 meV, respectively.For the same sphere without the shell, the FD is 240.2 μm.In our core-shell structures composed of concentric spheres and co-axial cylinders, the core due to its significantly larger size influences the focusing more than the shell does.Obviously, we do not expect to accurately calculate the characteristics of mesoscale objects using geometrical optics.However, it is interesting to note that the trend is in agreement with the FDTD simulation results.As the shell index decreases with the increase of the E f , the FD decreases.TJ formation is a complex interplay of the incoming, refracted, and diffracted waves inside and around the micro-object.Generally, smaller spheres and cylinders produce PJs closer to the distal interface, i.e. shorter FD, with smaller FWHM; also, higher refractive index of the core, would lead to shorter FD and smaller FWHM of the produced PJ 2 .In principle, these parameters can be further exploited to optimize the TJs.
In our previous work for 2D cases, we observed intensity enhancement factors in the range of ~ 9.5-22, here for the 3D cylinders the intensity enhancement factor is ~ 41-52.In addition to the refracted waves and diffracted waves at the periphery of the cylinder, the diffracted waves from the sides also contribute to the focal point's intensity in the 3D case. (2) Figure 3.The electric field distribution inside and around a BDS-coated dielectric sphere with R = 500 μm and d DS = 1.5 μm at T = 300 K for (a) E F = 10 meV (TE), (b) E F = 55 meV (TE), (c) E F = 10 meV (TM), and (d) E F = 55 meV (TM). (e) TJ's calculated FD (square symbols) and intensity (circle symbols) versus E F . (f) TJ's length (square symbol) and FWHM (circle symbol) versus E F .Hollow (solid) symbols correspond to TM (TE) polarization of the incident uniform beam.
Figure 5 shows how the optical properties of a TJ produced by a cube changes as the E F of the BDS coating changes.The distribution of the electric field of the TJs, for the TE-polarized incident wave, at E F = 10 meV and E F = 55 meV are shown in Fig. 5a,b, respectively.The quantitative properties of the TJs for the TM-polarized light source as a function of the E F are presented in Fig. 5c,d showing a similar trend as in the sphere and cylinder cases.However, the intensity in the cube-generated TJs is less than those produced by cylinders/spheres and the focal distance is longer due to the weaker focusing effect in the cube.The intensity enhancement of the TJ varied between 6.2 (6.4) and 8.6 (9.2) for TM (TE) case.The FD varies between 436.5 (458) μm and 397.5 (417) μm corresponding to 39 (41) μm tuning for TM (TE) case.Within the same Fermi energy, the FWHM with a tuning of 4.2 μm (5.6) reduces from 62.3 (56.4) μm to 58.1 (50.8) μm, while the length of TJ with a tuning of 49.4 (40) μm increased from 320 (342) μm to 369.4 (382) μm for TM (TE) case.The TE case provided more dynamic range for the tunability compared to the TM case for the intensity, focal distance, and the spot size.The TJ length showed more dynamic range for the TM case compared to the TE case.
In most of the cases, the focal distance and TJs are longer, more intense with smaller FWHM in TE cases compared to the TM cases.When comparing the findings for all 3D TJs, the highest tuning of the intensity is attained for spheres, while the tuning of the TJ's length is more noticeable for the cube.In all cases, as the E f increases, the intensity of the TJ increases.It can be in part explained by the fact that since the index of the BDS coating decreases as E f increases, the Fresnel losses at each interface decrease which increases the photon fluence at the focal point.The intensity of the TJ is highest for the sphere due to its extended curved boundary compared to the cylinder and cube that provides high photon fluence at the focal point.The TJ lengths are longer in the . (e) TJ's FD (square symbols) and intensity (circle symbols) versus E F . (f) TJ's calculated length (square symbols) and FWHM (circle symbols) versus E F .Hollow (solid) symbols correspond to TM (TE) polarization of the incident uniform beam.cubes (cylinders) than those in the cylinders (spheres).A quick convergence of the Poynting vector around the focal spot would produce more intense TJs with shorter length 95 .However, when the Poynting vector is less diverging and remains parallel near the focal spot, the TJ is more elongated with longer FD as in the cube case.
Figure 6 shows the TJ's electric field distribution (for E F = 55 meV) in the beam's eye-view, i.e. yz-plane encompassing the locus of maximum intensity.The TJ's FWHM along the y-axis is what we have reported so far.www.nature.com/scientificreports/The TJ's FWHM was calculated along the z-axis as well for these cases presented in Fig. 6.The FWHM values along y-axis (z-axis) are 31.2(25.7), 29.4 (69.7), and 58.1 (62.4) μm for the sphere, cylinder, and cube, respectively, under normal incidence and TM polarization of the uniform beam.The corresponding values for the TE case are 26.4 (32.8), 21.3 (49.3), and 50.8 (66) μm.It can be seen that the TJ distribution is not totally symmetric, a consequence of vectorial nature of the incident uniform beam.According to Richards and Wolf diffraction theory 96 , under high numerical aperture focus, the intensity profile of the focal spot would be elongated along the direction of the polarization.Here we see that effect in the produced TJs; for TE (TM) case elongation of the TJ can be seen along z-axis (y-axis).
Next, we analyzed the optical properties of BDS-coated dielectric objects under oblique incidence to verify the robustness of the produced TJs at various angles of the incident beam.We found that the intensity of the TJ was maximum at E F = 55 meV within the range of studied Fermi energies.Figure 7a-d illustrate the simulation results of electric field distribution in the xy-plane for source rotation angles of 5°, 25°, 40°, and 60° for the TE case.The corresponding results for the TM case are presented in Fig. 7e-h.When the source angle was altered, the TJ was deflected, and the intensity at the focused region was maintained over a large angular range, demonstrating the micro-object's ability to work as a wide-scanning TJ-focusing lens.When the angle of the incident uniform beam was varied, the numerical results for deflection angle of the TJ produced by the sphere is presented in Fig. 7i.The intensity as a function of the rotation angle of the incident uniform beam is shown in Fig. 7j.It should be noted that the uniform beam source was rotated along the x-axis which caused the asymmetric appearance of the TJs.It should be noted that the asymmetry is caused by the finite size of the incident uniform beam that leads to partial illumination of the object.Interestingly, it can be seen that the maximum intensity of the TJ occurs at ~ 5° source rotation angle, not necessarily at normal incidence.Introduction of an asymmetric optical path length along the optical axis due to the obliquity of the incident beam produces an asymmetric TJ as a result of difference between phase velocity and the interference of the waves inside the particle.
The simulation results for the cylinder case at E F = 55 meV are presented in Fig. 8 indicating how the TJ is deflected by changing the source rotation angle.Figure 8a-d  for input angles of 5°, 15°, 25°, and 40° for TE case.The corresponding results for the TM case are presented in Fig. 8e-h.Figure 8i shows the dependency of the deflected angle of the TJ versus the source rotation angle.The intensity as a function of the source rotation angle is shown Fig. 8j.The observed asymmetry in the TJ shape is due to the fact that the source was rotated along the x-axis.Similar to the sphere case, it can be seen that the maximum intensity of the TJ occurs at a non-trivial source rotation angle, here at ~ 15° source rotation angle.
The simulation results for the 3D BDS-coated cube at E F = 55 meV is shown in Fig. 9. Figure 9a-d show the electric field distribution for input angles of 5°, 10°, 20°, and 35° for the TE case.The corresponding results for the TM case are presented in Fig. 9e-h.Figure 9i,j show the relationship between the output angle and intensity as a function of the input angle.Similar to the sphere and cylinder cases, it can be seen that the maximum intensity of the TJ occurs at a morphology-dependent source rotation angle, here at ~ 10°.
It can be seen in Figs. 7, 8, 9 that the exit angle of the TJs in the sphere and cylinder cases follow a somewhat symmetric trend with respect to the incident angle due to the inherent circular symmetry of the producing micro-object.However, the asymmetry between θ in and θ out in the cube case is mostly stems from the lack of rotational symmetry of the cube compared to the sphere/cylinder.
A list of PJ characteristics produced using different objects at different wavelengths is provided in Table 1.These parameters should scale with the wavelength following Maxwell's equations.A distinct feature of the current work with the majority of the previous works is that the PJs are dynamically tunable.Depending on the application, operation wavelength, and choice of available materials, an ad hoc simulation must be performed to select the optimum PJ-producing geometry.
The proposed structures made of Cd 3 As 2 may be associated with certain considerations.For example, due to the high toxicity of Cd and As, extra caution is required in handling of these elements during the synthesis.Although various fabrication methods for Cd 3 As 2 have been reported, a high-yield cost-effective large-scale fabrication method has not been established yet, unlike the case for semiconductor chip fabrication.

Conclusions
Through FDTD numerical simulation, we demonstrated formation of tunable 3D photonic jets in terahertz spectral domain (terajets, TJs) by dielectric spheres, cubes, and cylinders coated with a BDS layer.The dependency of the TJ's main characteristics (FD, FWHM, intensity, and length) on the BDS's Fermi energy and incident uniform beam's polarization was explored for each case.Polarization of the illuminating light was found to play an important role on TJ characteristics.Simulation of the TJs under oblique incident demonstrated formation of asymmetric focused beams with the maximum intensity of the TJs at a non-trivial morphology-dependent source-angle.

Figure 1 .
Figure 1.Schematic of the simulation geometry.A uniform beam, propagating along the x-axis, impinges upon a sphere (radius R), a cylinder (radius R and height h), and a cube (length R) with refractive index n 1 coated with a BDS with thickness d DS and index n DS .(a) Terajet characteristics include focal distance (calculated from the proximal surface of the particle to the locus of the maximum intensity), beam FWHM in xy plane, TJ length (measured as the distance between distal to proximal 1/e intensity contour along the direction of the TJ, i.e. x-axis), and the maximum relative intensity.Geometry of oblique incident for a (b) sphere/cylinder and (c) cube.In the core-shell structure shown in panel (a), n c (R c ) and n s (R s ) represent the refractive indices (radii) of the core and shell, respectively.Drawings are not to scale.

Figure 2 .
Figure 2. The index of refraction of a BDS, mimicking Na 3 Bi or Cd 3 As 2 , versus E F at several THz frequencies at room temperature: (a) real and (b) imaginary part.

Figure 4 .
Figure 4.The distribution of the electric field inside and around a BDS-coated cylinder with R = 300 μm and d DS = 1.5 μm at T = 300 K for (a) E F = 10 meV (TE), (b) E F = 55 meV (TE), (c) E F = 10 meV (TM), and (d) E F = 55 meV (TM). (e) TJ's FD (square symbols) and intensity (circle symbols) versus E F . (f) TJ's calculated length (square symbols) and FWHM (circle symbols) versus E F .Hollow (solid) symbols correspond to TM (TE) polarization of the incident uniform beam.

Figure 7 .
Figure 7. (a-d) The electric field distributions of the TJ formed by a BDS-coated dielectric sphere with R = 500 μm and d DS = 1.5 μm at Fermi energy of E F = 55 meV for input angles of (a) 5°, (b) 25°, (c) 40° and (d) 60° for the TE case.(e-h) The TM counterparts of (a-d).The dashed circles show the boundary of the sphere.(i, j) The emergence angle and intensity of the TJ, respectively, as a function of the uniform beam rotation angle.

Figure 8 .
Figure 8. (a-d) The electric field distributions of the TJ formed by a BDS-coated dielectric cylinder at Fermi energy of E F = 55 meV for source rotation angles of (b) 5°, (c) 15°, (d) 25°, and (e) 40° for the TE case. (e-h) The TM counterparts of (a-d).The dashed circles represent the boundary of the cylinder.(i, j) The emergence angle and intensity of the TJ, respectively, as a function of the uniform beam rotation angle.

Figure 9 .
Figure 9. (a-d) The electric field distributions of the TJ formed by a BDS-coated dielectric cube with R = 300 μm and d DS = 1.5 μm at Fermi energy of E F = 55 meV for source rotation angles of (a) 5°, (b) 10°, (c) 20°, and (d) 35° for the TE case.(e-h) The TM counterparts of (a-d).The dashed squares represent the boundary of the cube.(i, j) The output angle and intensity of the TJ, respectively, as a function of the uniform beam rotation angle.

Table 1 .
Characteristics of photonic jets produced by objects with various shapes and sizes.DS Dirac semimetal, LC liquid crystal, N/R Not reported, n c core index, n s shell index, n cl cladding index.